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The efficient placement of signal templates in source-parameter space is a crucial requisite for 
exhaustive matched-filtering searches of modeled gravitational-wave sources, as well as other searches 
based on more general detection statistics. Unfortunately, the current placement algorithms based 
on regular parameter-space meshes are difficult to generalize beyond simple signal models with few 
parameters. Various authors have suggested that a general, flexible, yet efficient alternative can 
be found in randomized placement strategies such as random placement and stochastic placement, 
which enhances random placement by selectively rejecting templates that are too close to others. 
In this article we explore several theoretical and practical issues in randomized placement: the 
size and performance of the resulting template banks; the very general, purely geometric effects 
of parameter-space boundaries; the use of quasirandom (self-avoiding) number sequences; most 
important, the implementation of these algorithms in curved signal manifolds with and without the 
use of a Riemannian signal metric, which may be difficult to obtain. Specifically, we show how the 
metric can be replaced with a discrete triangulation-based representation of local geometry. We 
argue that the broad class of randomized placement algorithms offers a promising answer to many 
search problems, but that the specific choice of a scheme and its implementation details will still 
need to be fine-tuned separately for each problem. 

PACS numbers: 04.30.Db, 04.25.Nx, 04.80.Nn, 95.55.Ym 

I. INTRODUCTION 

An international network of ground-based laser- interferometric gravitational- wave (GW) detectors [IH3]; which 
operate in the high-frequency band between 10 and 10 3 Hz, has by now completed several rounds of science data 
collection and analysis. Space-based detectors such as LISA [4] will extend our reach to frequencies as low as 10~ 5 
Hz, and to sources as far as redshifts of 20. The first direct detections of GW sources will inaugurate the era of 
GW astronomy, a novel mode of inquiry of the Universe that will provide unique information about the properties of 
relativistic objects such as black holes and neutron stars, and allow precise tests of general relativity's yet unproven 
predictions. 

Because GW signals come immersed in noisy detector data, the statement that a GW source was detected is by its 
nature probabilistic, and can be quantified with frequentist false-alarm and false-dismissal probabilities jS], or with 
a Bayesian odds ratio [6] . In an ideal world where the statistical properties of detector noise are known exactly and 
detector calibration is perfect, detection corresponds to a simple mathematical statement about the data [7], although 
actually evaluating that statement may require considerable numerical computation. In the real world where noise is 
non-Gaussian, nonstationary, and poorly characterized, and where calibration is challenging, detection corresponds 
to a more complicated, heuristic process, which involves expert judgment (such as the choice of data-quality vetoes 
and multidetector coincidence windows), and which is embodied in complex software pipelines [SJ. The validity 
and efficiency of the process are evaluated empirically, by measuring the resulting background of spurious detection 
candidates, and by injecting simulated GW signals and verifying their recovery. Indeed, the day-to-day work of GW 
analysts is very much dominated by the development and validation of such pipelines. 

Searches for GW signals of known shape from well-understood sources, such as the archetypical gravitational 
two-body problem of compact-binary inspiral, merger, and ringdown, occupy a privileged spot in the GW data- 
analysis landscape because they are likely to yield the first robust detections, although searches for unmodeled GW 
bursts in coincidence with EM counterparts j9] are also a contender. Furthermore, modeled-signal searches find 
a straightforward mathematical realization in the simple yet powerful technique of matched filtering [5]: roughly 
speaking, the detector output is time-correlated with a set of theoretical waveforms (templates), and a peak in one 
of the correlation time series indicates that a GW signal of that shape may be present in the data at the time of the 
peak. Higher peaks correspond to higher detection confidence, because stronger signals have smaller probability of 
being simulated by random features in the noise. 

The template placement problem arises because we are interested in searching for GW sources with a range of phys- 
ical parameters (such as the component masses in a binary), and because different parameters can yield considerably 
different waveforms. Therefore we must choose a finite number of templates to correlate with the detector data; these 
templates must span the region of interest in source parameter space; they must be placed "densely" enough to avoid 
reducing detection confidence (the height of the correlation peaks) by matching a true GW signal with too different 
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a template; yet they must be spaced as "sparsely" as possible to limit their total number and therefore the compu- 
tational cost of the search. The notion of distance implied here by "densely" and "sparsely" is made quantitatively 
precise by defining a noise-weighted distance A [Eq. <t3J) in Sec. O in the linear space M. Nbi * of all possible detector 
outputs. The distance A between two GW signals is directly related to the reduction in detection confidence incurred 
by using one to match the other, and the noise weighting encodes the different sensitivity of the detector at different 
frequencies. Thus we wish to place templates uniformly from the "point of view" of the detector. 

Now, a family of signal templates {h(Xj)} parametrized smoothly by the source parameters Xj realizes an 2V Bma u- 
dimensional topological manifold H embedded in R-^* . Thus the placement problem consists in selecting a finite and 
discrete bank {h(X^)} C % such that the (K Arbis -)distances of all possible incoming GW signals (i.e., all the points in 
H) from the nearest bank template h(X^) satisfy certain properties. For instance, we may constrain the average and 
variance of this distance, or we may require that the distance may be no larger than a maximum value. In the latter 
case, as emphasized by Prix [10] . template placement becomes an instance of the sphere-covering problem: given the 
manifold H, we wish to cover it with the smallest possible number of spheres of fixed (A) radius, where the sphere 
centers sit on % and correspond to the {h(X^)} in the bank. Qualitatively speaking, we want the spheres to cover 
the manifold completely, so that we miss no signal; at the same time we wish to reduce their overlap, so that we not 
waste computational resources on testing templates that are almost duplicates. 

The distance A endows H with a local Riemannian structure, and therefore an iV sma n-dimensional metric tensor 
9jk [11H13] . Template placement can then be reformulated as the problem of covering % with spheres of fixed geodetic 
radius. If T~L is curved, the geodetic distance will consistently overestimate the A distance, so the reformulation is 
only faithful when curvature is negligible over the typical sphere radius. GW analysts have generally been happy to 
accept this approximation and thus work in the familiar framework of differential geometry; so shall we. 

If % has no intrinsic curvature, and if a suitable change of coordinates can be found that maps H to Euclidean 
space, 1 therefore mapping g^ to the identity matrix, we can rely on a wealth of theoretical results about the covering 
properties of periodic lattices, and place sphere centers (i.e., bank templates) at the lattice vertices: for instance, on 
a square [17] or hexagonal [HI [19] lattice on the Euclidean plane, or on their generalization in higher-dimensional 
Euclidean space [ID] . If however no such change of coordinates can be found, perhaps because the metric itself is 
difficult to obtain or unreliable, or if % has significant curvature, it becomes very hard to place templates along any 
regular structure. In the case of two source parameters, Arnaud et al. [20] and Beauville et al. [H] describe algorithms 
that create curvature-conforming tilings by propagating a locally hexagonal lattice while adjusting for the change in 
the equal-distance contours. Such an approach has not been attempted in higher-dimensional spaces, where it appears 
very cumbersome. 

A possible alternative is placing templates at the vertices of suitably "grown" unstructured meshes, such as adap- 



tively refined triangulations. We investigate such a strategy in a separate article [22], and briefly in Sec. VII below. 
Here we concentrate instead on the randomized sampling strategy originally proposed by Allen [23j and later studied 
by Babak [24] , Messenger et al. [25], and Harry et al. [2B]. In its simplest incarnation (pure random sampling), the 
idea is to draw a set of points that are (pseudo-)randomly distributed across W. Because of its very randomness, 
such a covering can never be guaranteed to cover 100% of T-L\ however, as shown by Messenger and colleagues [25], 
it can achieve very high covering fractions with total numbers of points that (in Euclidean space) are competitive 
with the most efficient relaxed lattice coverings (i.e., lattices rescaled to cover a fraction < 100% of space). To cover 
curved manifolds, it is sufficient to know the metric determinant, and adjust the local probability of random draws 
accordingly. (Recently, Rover has suggested that the Bayesian prior distribution for the source parameters can be used 
to adjust the local template density in order to minimize the expected number of draws until a matching template is 
found [27].) 

In this paper we extend the results of Ref. [25]: we investigate the finer details of random coverings, such as the 
variance of the covering fraction and the effects of boundaries, and we provide further evidence that uniform coverings 
of curved manifolds can be obtained in general conditions, even without knowledge of the metric, but using only the 
A distance. Specifically: 



We confirm that random coverings are more efficient than could be naively expected, and we provide an exact 
formula for the variance of the covering fraction. 

We extend Messenger et al.'s in-the-bulk results by considering the effects of manifold boundaries in reducing 
the covering fraction: we find that these effects become overwhelming, and must be treated carefully, as the 
manifold dimension increases. For example, neglecting boundary effects, 11,000 random spheres of radius 0.3 



1 For instance, approximate but accurate mappings can be found for the 2.5PN restricted inspiral signals from binaries of spinless compact 
bodies 1141 1151 and for the continuous GW signals from rapidly rotating neutron stars |16| . 
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cover the 8-dimensional unit cube with 95% covering fraction; including boundary effects, 92,000 spheres are 
needed. 

• We show that quasirandom numbers (self-avoiding deterministic sequences often used to reduce variance in 
Monte Carlo integration |28j ) can be used instead of pseudorandom numbers to improve the efficiency of random 
coverings, although such a substitution is most effective on quasiflat manifolds. 

• We investigate stochastic placement strategies [Ml where the random distribution is tempered by considering 
the distances between selected points (e.g., rejecting draws that are too close to another point in the bank): we 
quantify the efficiency of these coverings, and examine their application to curved manifolds without the use of 
the metric. 

• Last, we propose an alternative metricless placement strategy, based on triangulating the curved manifold to 
represent its local geometry as the set of finite distances between neighboring triangulation points. 

Overall, we find that even flexible placement strategies such as random, stochastic, and unstructured-mesh coverings 
cannot be made into recipes of general applicability, but need to be tailored to the circumstances of different template 
families and searches. Much as the simple mathematical formalism of matched-filtering GW searches has made way 
to complex processes and pipelines, so must template placement ultimately evolve from a problem of mathematics to 
one of engineering. On the other hand, we expect that the general methods and results discussed in this paper would 
apply also to the larger class of GW searches that use more general detection statistics than the matched-filtering 
SNR (j2j, but that still require the placement of (generalized) templates according to notion of distance. Examples 
include stack-slide , Hough-transform [3D], J"-statistic [31], and global-correlation [TB] searches for periodic GWs. 

We note also that not all matched-filtering searches rely on template banks. The broad family of Monte Carlo 
methods, which include Markov Chain Monte Carlo, genetic algorithms, and nested sampling (see [32] for a LISA- 
centric review) , conceptualize searches as the randomized exploration of template-signal correlation across parameter 



space, rather than its evaluation at a set of predetermined points. However, as we argue in Sec. VIII the systematic 
exploration of parameter space enabled by random and stochastic banks, as well as the scaling properties derived here 
and in Refs. [25] [26], may find applications in designing or guiding Monte Carlo searches. 

This paper is organized as follows. In Sec. [H we recall key formulas from the theory of matched filtering, and 
emphasize their geometric interpretation. In Sec. III[ we introduce the sphere-covering problem. In Sec. IV we study 
the properties of random coverings. In Sec. [V] we examine the effects of boundaries. In Sec. |VI[ we quantify the 
improved efficiency possible with quasirandom sequences and stochastic placement. In Sec. |VII[ we explore the use 
of unstructured meshes for metricless placement. In Sec. |VIII[ we summarize our conclusions, and discuss future 
directions of investigation. In the Appendix, we describe the numerical procedures used throughout this work. 



II. KEY MATCHED-FILTERING FORMULAS 



Matched-filtering searches for GW signals are based on the systematic comparison of the measured detector output 
s with a bank of theoretical signal templates {h(\^ k ')}. Here we develop only the formulas that we need later, but a 
few useful starting points in the vast literature on this subject are Refs. [51 [71 133TI37| . 

Inner product. A crucial mathematical construct is the symmetric inner product ((71,32) between two real signals 
(detector outputs over a fixed time period) 171 (t) and 32 (t) ■ This product is essentially the cross-correlation between 
gi(t) and 32 (t), weighted to emphasize the frequencies where detector sensitivity is better. We follow Cutler and 
Flanagan's conventions [36 j and define 

(31,52) =2/ o~7T7h J ~ / £T77\ d /' C 1 ) 



Sn(\f\) ' Jo S n (f) 

where the tildes denote Fourier transforms, the stars complex conjugates, and S n (f) the one-sided power spectral 
density of detector noise. We can then define the signal-to-noise ratio (SNR) p of detector output s after filtering by 
template h(X^) = as 

o(h^) = {s > kik)) = {s > k(k)) (2) 
Hy ' rms(n,/i«) y/(hV°),h( k ))' 

where n denotes detector noise, and ims{n,h^} represents the rms average of that inner product over all possible 
noise realizations (see, e.g., [5]). It is convenient to think of Eq. [2] as the SNR of s after filtering by the normalized 
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template = hffl j-\J (h( k \ ft,( fe )). If s consists solely of Gaussian and stationary detector noise, then p is a normal 
random variable with mean zero and variance one. If instead a signal equal to A is present in s in addition to noise, 
then p is a normal random variable with mean A and again variance one. We claim a detection of a signal with source 
parameters A' fc J whenever p{h^ k ') is greater than a chosen threshold p* . Such a detection scheme is Neyman-Pearson 
optimal in the sense that it maximizes the probability of correct detection for a given probability of false detection. 

Distance. Henceforth, we shall restrict our consideration to normalized signals, and drop the caret that we just 
used to denote them. Given two such signals, we can use the inner product to define a distance between them, 



g 2 ) = Vwi - .92,. 9i -52)/2 = V 1 - (siwffa)- (3) 

This distance is a proper metric, since it is semidefinite positive (and vanishes only for g\ — g 2 ), symmetric, and it 
satisfies the triangle inequality, 

A(g 1 ,g 2 )<A(g 1 ,g 3 ) + A(g 3 ,g 2 ). (4) 

As mentioned in the introduction, we think of signals as points in the metric space M. Nbis of all possible detector 
outputs; A^big is large but finite, because the response of the detector and the frequency content of prospective GW 
signals are band-limited, and thus can be represented by a finite number of samples (or equivalently Fourier coefficients, 
or coefficients over some other basis). Furthermore, a signal family {h(X)} (e.g., the GW signals from neutron-star 
binaries over a given range of masses) that is smoothly parametrized by the source parameters Xj realizes an iV sma ii- 
dimensional topological manifold embedded in R bl «. The inner product endows this manifold with a Riemannian 
structure, and the distance between nearby points on it is then approximated by a metric tensor g kl , 



1 d 2 (/i(A 7 ),/i(A')) 

A 2 (h(Xj), h(Xj + dX j ))m - 



2 dX' k dX[ 



dXkdXi = gki(Xj)dXkdXi. (5) 



A, 



One criterion to choose the bank templates {h(X^)} from the continuous family {h(X)} is to require that for every 
signal in {/i(A)}, the nearest bank template be at most at a distance A max . This ensures that the SNR is reduced at 
worst 2 by a factor MM = 1 — A^ ax , known as the minimum match. 3 With this criterion, template placement can be 
seen as an instance of the covering problem, the subject of the rest of this paper. 

Extrinsic parameters. It is sometimes possible to factor out certain source parameters (known as extrinsic) from 
the template placement problem: this happens when the SNR can be maximized analytically or numerically over the 
appropriate range of the extrinsic parameters, so it is not necessary to explicitly place templates over a discrete set of 
their values. Two examples are the time of arrival and initial phase of binary inspiral signals. Extrinsic parameters 
are usually welcome, because they reduce the computational cost of searches; on the other hand, they complicate the 
geometric interpretation of template placement, because the distance between two signals is normally a function of 
all source parameters, so we need a prescription to choose representative values of the extrinsic parameters. 

If we denote the extrinsic parameters collectively as </>, and the intrinsic (nonextrinsic) parameters as 9, the logic 
of matched-filtering detection would suggest using the minmax prescription 



A{h{0i),h{9 2 )) = Jl- minmax {h{e 1 ,cj> 1 ),h(e 2 ,h)) ■ (6) 



if we think of h(0i,(f>i) as the GW signal to be detected, and of h(0i, <j>±) as the nearest template, then max^ 2 represents 
the "automatic" maximization of SNR over the extrinsic parameters, while min^ would conservatively ensure that 
the minimum match is achieved even in the least favorable case. Unfortunately, Eq. Q does not yield a true distance: 
it is not symmetric, and it does not satisfy the triangular inequality. 

An alternative is to choose the extrinsic parameters as arbitrary smooth functions of the intrinsic parameters, 

A(h(e 1 ),h{e 2 )) = Vi-(/i(8i^(8i)),/i(fl2,^(«2))), (7) 

and verify that it yields an acceptable distribution of (^-maximized SNRs over all <pi . (For instance, we may choose 
4>(0) to locally maximize the determinant of the full metric gjk-) Throughout the rest of this paper, we shall assume 
that a sensible distance along the lines of Q can be defined whenever extrinsic parameters are present. 



2 Because SNR is inversely proportional to luminosity distance, the horizon distance to which a given source can be detected is reduced 
at worst by a factor MM. It follows that if sources are distributed uniformly around the detector, the rate of event detection will be 
reduced by a factor not worse than MM 3 . 

3 However, at least in the case of binary inspiral signals, interpolation across a sparse cubic-lattice bank can be used to retrieve a fraction 
of SNR larger than the bank's nominal MM 15, 38 139] . The gains possible with this technique are not reflected in our analysis, which 
is concerned with the general covering problem rather than with specific cases. 
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TABLE I: Volume of the d-dimensional unit-radius sphere. Interestingly, Vd is maximum for d = 5. 



III. THE COVERING PROBLEM: PERIODIC LATTICES 



Prix |10j discusses template placement as an instance of the sphere-covering problem. See that paper and Ref. [25 
for a general introduction. Here we focus on essential formulas and new insights. We adopt the notation of Ref. |25j . 

Thickness. The problem of covering d-dimensional Euclidean space E d with a regular periodic lattice is well known 
in the mathematical literature [40] : it consists in finding the lattice with the smallest number of points per unit volume 
such that, if a <i-dimensional sphere of fixed radius r is placed with its center at each lattice point, the spheres cover 
E d completely. A typical quantity used to evaluate the efficiency of such a covering is the thickness 0, defined (for a 
subvolume V m of E d ) as 

N Vd r d 

6 = ^7^, (8) 

where N is the number of lattice points (and therefore spheres), Vd the volume of the d-dimensional unit-radius sphere, 
and r the required sphere radius. The unit-sphere volume V c i can be expressed in terms of the gamma function as 
V d = ir d/2 /T{d/2 + 1); its values for d up to 20 are given in Table Q 

The thickness is a pure number: it is the ratio between the total volume of the spheres used in the covering and the 
total covered volume. The best possible case is a covering with 9 = 1, where the spheres do not intersect, and each 
point of E d is inside one and only one sphere; of course it is impossible to reach this limit for d > 1. Mathematicians 
have worked out theoretical limits for the thinnest possible coverings, such as the CFR bound [?D] 

^ d d . . 

The most efficient known lattice coverings for d < 20 are listed with their 8 in columns 2 and 3 of Table [ITJ Note that 
the authors of Ref. |25j usually report values for the normalized thickness 9 = Q/Vd rather than for 0, as we do. 

Partial covering. A useful generalization of the problem is to allow a fraction of space to remain uncovered, and 
define an effective thickness that depends not only on the choice of lattice, but also on the percentage of the volume 
that it covers. In the language of template banks, this means that the maximum-distance criterion will not be 
achieved for a fraction (usually small) of the possible incoming GW signals. As noticed in Ref. [25], the effective 
thickness of partial regular coverings (even for rather large covering fractions) improves dramatically compared to 
complete coverings, especially as d increases. This means that in a complete covering most of the lattice density is 
needed to reach the "hardest" few percent of the volume; by contrast, even a simple 1 d (simple hypercubic) lattice 
can efficiently cover 95% or 99% of space. See columns 4-11 of Table [Tl] for the partial-covering thicknesses of Z d 
and (Ac?)* lattices. (The latter are the d-dimensional generalization of body-centered cubic lattices, whereas the Ad 
lattices are the generalization of face-centered cubic lattices. Here the star denotes a reciprocal lattice. A2 is the 
hexagonal lattice, and is its own reciprocal. See Ref. [JD] for details.) 

In this context, a better index of covering efficiency is not the complete-covering thickness @ioo%, but the variation 
of as a function of the covering fraction. The left panel of Fig. [ljshows this function for (Ad)* lattices with d = 2-8. 
The right panel shows the isotropic dilation factor that must be applied to a complete covering to achieve a lower 
covering fraction with the same thickness. Thus we could say (for instance) that an (A8)* covering is ~ 2.3 times 
more efficient at an 80% covering fraction than at 100%, where the factor 2.3 is the ratio of the thicknesses, which is 
also equal to the eighth power of the dilation factor. 

Curved manifolds. So far we have considered the covering problem for Euclidean space, and equivalently for flat 
manifolds H that can be mapped to K d by a coordinate change. What about curved manifolds? For these, the metric 
gki is a function of position on the manifold, and the efficiency of lattice coverings that are periodic with respect to 
the coordinates drops considerably. This is because of three different effects, illustrated in Fig. [2] the variation of 
the determinant of the metric, the variation in the ratio of its eigenvalues, and the variation in the orientation of 
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20.37 
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TABLE II: A comparison of thickness as function of d for the most efficient known complete coverings (column 2), for partial 
and complete (Ad)* coverings (columns 4-7), for partial and complete Z d coverings (columns 8-11), and for random coverings 
(columns 12-15). In this table, thicknesses are computed in the limit V m — > 00, neglecting any boundary effects. See also |41ll42j . 
All values were computed using the technique discussed in the Appendix, except for those in italic, which are taken from Ref. 
[25] . and reported here for completeness. 

its eigenvectors. This distinction is somewhat academic at this point, but it will acquire a deeper meaning in later 
sections. 

Let us first consider the effects of metric-determinant variation. We may refer to the right panel of Fig. [T] to gain a 
quantitative understanding in the case of (Ad)* lattices. Suppose that one such lattice was arranged to provide 100% 
covering around the point where \gu\ is minimum: around a different point A i? the lattice will be in effect dilated 
isotropically by a factor {\gki{\)\/\gki{\)\) 1 ^ 2d \ the corresponding reduction in the covering fraction can be read off 
from Fig. [T] Changes of a few tens % in \gu\ are already very damaging. 

Let us now consider the effects of variation in the ratio of metric eigenvalues. This can be interpreted as a 
local dilation and contraction of space along orthogonal axes, as displayed in the left panel of Fig. [3] for an (A2)* 
(hexagonal) covering. The covering becomes insufficient along one direction, and inefficient along the other. The 
right panel of Fig. [3] shows how the dilation-contraction affects the covering- fraction-thickness curve. The worsening 
covering performance can always be seen from two opposite viewpoints: either as a loss of efficiency (i.e., an increase 
in thickness) for the same covering fraction, or a decrease in covering fraction for the same thickness. 

In fact, the thickness of a periodic lattice can be made arbitrarily large even while keeping \gu\ unchanged. The 
reason is that the points of periodic lattices lie on hyperplanes, since they can all be obtained as linear combinations 
with integer coefficients of d independent vectors [ID]. Consider now the transformation given by a simultaneous 
contraction inside the hyperplanes defined by d — 1 of those vectors, and a dilation in the orthogonal direction. 
(Arrange the dilation-contraction factors to cancel out in \g^i | .) Much like what happened in Fig.|TJ the transformation 
increases the sphere radius needed to cover the space between hyperplanes, but also the superposition of the spheres 
along directions contained in the hyperplanes. 
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Covering fraction [%] Covering fraction [%] 

FIG. 1: Efficiency of partial coverings. Left: thickness of (Ad)* coverings as a function of the covering fraction, for d = 2-8. 
Right: isotropic dilation factor between a complete covering and partial covering of the same thickness, for the covering fraction 
given on the horizontal axis. 



Position-dependent Position-dependent Position-dependent 

metric determinant metric eigenvalues metric eigenvectors 




FIG. 2: Three ways in which manifold curvature disturbs a lattice covering laid out periodically along manifold coordinates. 
Left: variation in the metric's determinant causes spheres with the same geodesic-distance radius to have different coordinate 
radii. Center: variation in the ratio of the metric's eigenvalues maps geodesic-distance spheres into coordinate ellipses. Right: 
variation in the metric's eigenvectors changes the orientation of coordinate ellipses. Because of all three effects, the covering 
fraction comes to depend on the position in the manifold. 



IV. THE COVERING PROBLEM: RANDOM COVERINGS 



In Euclidean space, a random covering is obtained by choosing TV points randomly and uniformly distributed 
across the region to be covered. Such an arrangement can never be guaranteed to cover 100% of points in all of its 
realizations; random coverings however generate very efficient partial coverings of E d — indeed, for sufficiently high d, 
the most efficient partial coverings [25] . 

Effective thickness. We can use probabilistic reasoning to characterize for a random covering. Assume that the 
probability of choosing any point in E d is the same; given a (generic) reference point P, the probability that a sphere 
of radius r with randomly chosen center will contain P is 

Pin = V d r d /V m < 1, (10) 

the ratio of the volume of a sphere to the manifold volume to be covered. The probability that N randomly centered 
spheres will not contain P is then 

Vn = (1 - V d r d /V m ) N < 1, (11) 

so the probability that at least one of the N spheres covers P is pc = 1 — Pn- Taking the logarithm, if the spheres 
are small enough compared to V m (more precisely, for V m and N — > oo with constant V m /N), 



log(l-pc) = AHog(l - V d r d /V m ) ~ -NV d r d /V m = -9 



(12) 
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Covering fraction [%] 

FIG. 3: Left: effect of dilating space along the x direction and contracting along y, for a d — 2 hexagonal covering. The 
dilation-contraction factor is 2, so the determinant of the transformation is 1. The spheres centered on the lattice points are 
now too far apart along the x direction, too close along y. The new lattice realizes a covering with the larger radius (i.e., 
maximum template distance) of the dashed circle. Notice that the Voronoi cells (i.e., the sets of points closer to a given lattice 
point than to any other lattice point) of the stretched lattice are not the stretched Voronoi cells of the original lattice. Right: 
reduction in covering fraction (for a fixed thickness) or increase in thickness (for a fixed covering fraction) for the dilation- 
contraction of the left panel, shown here also for a d = 2 cubic covering. As we shall see in Sec. |IV| random coverings are not 
affected. 



[according to definition |(8j|]. Thus the effective thickness is a function of the probability that a generic point is covered, 

% a a dom - lo s (r^) ■ as) 

As shown in the rightmost "Random" section of Tablc|nj 0_^° m = oo, %l% ora = 4.6, 6^ om = 3.0, and e^ om = 
2.3, independently of the dimension d. These numbers assume the limit r — > 0, not just because of the approximation 



in Eq. (12 1, but also to cancel out boundary effects. As we shall see below, these become critical for large d, or when 
one or more dimensions have extents comparable to r. Strictly speaking, ©g^ 0111 is a function of pc (the probability 
that a generic point is covered), and not of the covering fraction C, which is a random variable that depends on the 
particular realization of the covering. However, we shall see shortly that the average of C over realizations is just pc, 
while its variance vanishes in the limit of N — > oo (and therefore of r — > 0) . 



Equation (13 1 yields an estimate of the number of random points needed to achieve a given pc- For instance, for 
Pc = 95%, 

» = ear (£) - * (r^k) (M * 3 (£) • <") 

Let us also establish (for future use) an expression for the probability that at least one of N randomly chosen points 
falls inside a region of volume Vr. Following the same line of reasoning as above, for N ^> 1 we find 

rm.N i (, © Vr \ N _^ _____ ( ^ V R 



= 1 " ^ " vr/V} = 1 -{ 1 -nv^) * 1 - exp lr e v# ) ■ (15) 

Average and variance of the covering fraction. Following Messenger and colleagues [25] . let us now confirm that 
the average covering fraction is pc\ extending their derivation, we shall also obtain an exact expression for the variance 
of the covering fraction. 

We define the characteristic function f(P) to be one for points that are covered in a realization of the covering, and 
zero for points that are left uncovered. The covering fraction is then C = (1/V m ) J f(P)dV m (where P is parametrized 
by Aj, and dV m = \/\gki(K)\dXi) and its expectation value is 

E \ C \ = TF~ I E[f{P)]dV m = -J- / Pc dV m = pc- (16) 



V / LJ v /j v 



The computation of the variance is more involved: 
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FIG. 4: Left: typical geometry that enters the computation of the variance of the covering fraction and of boundary effects. 
The basic ingredient is the volume of the slice of a sphere cut by parallel hyperplanes. Right: boundary effects. The probability 
that a generic point is covered scales with the "available" manifold volume within a distance smaller than the covering radius. 



E[C 2 } - E[C] 2 = ^jj E[f(P)f(P')]dV m dV^ - P 2 C = ±J (E[f(P )f(P)] - p%) dV m . (17) 

If points Pq and P are separated by a distance larger than 2r (twice the covering radius), no sphere can cover both, 
so the probability that P is covered and the probability that P is covered are uncorrelated, E[f(P )f(P)] = p 2 c . 
Switching to a set of Euclidean, spherical coordinates centered in Pq, we then have 



E[C 2 } -E[C} 2 = ^[ p »-i (E[f(0)f(p)] - p%) dp, (18) 

Vm .1 p<2r 



where Sd is the surface of the <i-dimensional unit sphere. We define Pboth(p) = E[f(0)f(p)] to be the probability 
that two points separated by distance p are both covered. Looking at the left panel of Fig. |4j we conclude that 
Pboth = 1 — (1 — P2Z 2 ){^ ~ Pz)> where P2Z 2 is the probability of having at least one covering point in either of the 
regions labeled Zi , while p% is the probability of having at least one point in both the regions labeled Z\ . Now, the 
volumes of Z\ and Z2 can be expressed in terms of the volume covered by the portion of a d-dimensional sphere (of 
radius r, the covering radius) that is cut by parallel hyperplanes at (signed) distances hi and ft 2 from the center (see 
again the left panel of Fig. Eh, 

/>arcsin(/i2 / f~) 

V d l (hx,h 2 ,r) = V d -ir d / cos d (9)d8; (19) 

J arcsin(/ii fr) 

we then have 

VzM = V d Hp/2,r,r), V Zl (p) = r d V d -2V Za (p). (20) 



Using Eq. ( 15 ) to express P2Z 2 an d p Zl , after some manipulation we find 

f-n/2 



Pboth(p) = l-2cxp(-e) + cxp (-26 + 29^ f cos d {6)d0 j , (21) 

V *d Jarcsin(p/2r) I 

and finally 



E[C 2 } - E[C] 2 = ^ f y d - x [ Photh (ry) - p 2 c ] dy, ( 22 ) 







where we have used the fact that Q/N = Vdr d /V m and that Sd — dVd- Note that Eq. (22 1 does not depend on 
the covering radius r. Table |III| gives a few numerical values for this expression, which can be used instead of the 
"guesstimate" of Ref. E[C' 2 ] — E[C] 2 ~ PcO- ~ Pc)/(2dN). The two expressions share the initial decrease of 
variance with increasing d, and the dominant 1/N scaling. 

Maximum distance. Another obvious indicator of the performance of random coverings is the maximum distance 
r w of a point in the covered manifold region from the nearest covering point. This quantity coincides with the radius of 
the largest empty sphere that can be fit among the covering points, and again it is a random variable that depends on 
the particular covering realization. The authors of Ref. 25j make another guesstimate for the probability distribution 
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d 


50% 


80% 


90% 


95% 


99% 


d 


50% 


80% 


90% 


95% 


99% 


1 





1534 





1913 


0.1340 


0.0800 


0.0189 


11 





1212 





1059 


0.0548 


0.0234 


0.0023 


2 





1421 





1578 


0.1003 


0.0540 


0.0099 


12 





1210 





1053 


0.0543 


0.0232 


0.0022 


3 





1351 





1389 


0.0826 


0.0413 


0.0062 


13 





1207 





1049 


0.0540 


0.0230 


0.0022 


4 





1306 





1274 


0.0723 


0.0343 


0.0045 


14 





1206 





1046 


0.0537 


0.0228 


0.0022 


5 





1276 





1201 


0.0661 


0.0303 


0.0036 


15 





1205 





1043 


0.0535 


0.0227 


0.0022 


6 





1255 





1152 


0.0620 


0.0278 


0.0030 


16 





1204 





1041 


0.0534 


0.0227 


0.0022 


7 





1240 





1119 


0.0594 


0.0261 


0.0027 


17 





1203 





1040 


0.0533 


0.0226 


0.0021 


8 





1229 





1096 


0.0576 


0.0250 


0.0025 


18 





1203 





1039 


0.0532 


0.0226 


0.0021 


9 





1222 





1079 


0.0563 


0.0243 


0.0024 


19 





1202 





1038 


0.0532 


0.0225 


0.0021 


10 





1216 





1068 


0.0554 


0.0238 


0.0023 


20 





1202 





1038 


0.0531 


0.0225 


0.0021 



TABLE III: Variance of the covering fraction [Eq. (22 1], before scaling by 1/N. Clearly the covering fraction converges very 
quickly to pc for even moderate N. The dependence on d is mild. 



of r w : in a given realization of a random bank in d dimension, with N templates and a covering radius r the probability 
that maximum distance is smaller than f would be 



p(r w < f) ~ [1 - exp(9(f/r) d )] 



2dN 



(23) 



The guesstimate fits reasonably the results of low-dimensional numerical experiments with random coverings at 90% 
covering percentage. Anyway, since this probability distribution is based on the same approximation used for the 
guesstimate of the variance of the covering fraction, it must be used cum grano salis: it may not provide, as claimed, 
10%-accurate predictions at higher d, or with different covering fractions. Still, the Monte Carlo simulations presented 
in Ref. |25j do indicate that the mean and variance of r w decrease with increasing d, but it is not clear whether they 
tend asymptotically to zero (a complete covering in every realization). 

Curved manifolds. If we leave Euclidean space, the proper definition of "random" is that the probability for each 
point in the covering to lie within a manifold subregion should be proportional to the volume of that subregion. The 
derivation that led to Eqs. (13) and (14) can then be reproduced verbatim, with the caveat that r must be small 
enough compared to the curvature scale for Vd,r d to be a good approximation to the volume of a sphere of radius r 
(it is hard to do without this assumption, which we shall maintain to be valid in the rest of this paper). 

In coordinate space, random manifold points are distributed with local density proportional to v/[gfci|. To achieve 
this in practice, we may for instance use a rejection method whereby points A; are chosen uniformly in coordinate 
space, and then accepted with probability \f\gu (Aj) \/{vaax\ i \J\gkl{\)\)- The computational cost of this procedure 
will be dominated by the estimation of the total manifold volume through the integration of and by its 

local evaluation. Thanks to the modulation of (coordinatewise) density, and to the fact that random coverings are 
distributed along no preferred hyperplanes, and with no preferred lengthscales, random coverings on curved manifolds 
avoid all three problems illustrated in Fig. [2j and achieve the same thickness (for the same covering fraction) as in 
Euclidean space. 



BOUNDARY EFFECTS FOR RANDOM COVERINGS 



The analysis of Messenger and colleagues [25 , as well as our discussion so far, have dealt with the thickness 
of coverings in the bulk, and neglected boundary effects. Unfortunately, these can be very important, and can be 
difficult to correct in practice. For random coverings, it is relatively easy to understand their origin and estimate 
their magnitude. As illustrated for points A and C in the right panel of Fig. |4j at distances < r from the boundary, 
the probability that a generic point is covered is smaller than farther away in the bulk; this is because the probability 
scales with the manifold volume contained within a radius r from the point. 

Covering probability. To obtain the covering probability for the points on a boundary, we recall Eq. ( 15 ) for the 
probability of having at least one covering point in a region of volume Vr, and obtain pbound,d-i = 1 — exp(8/2) = 1 — 
(1— pc) 1 / 2 for the center of a half sphere on a (d—1) -dimensional boundary, pbound,d-2 = 1— exp(9/4) = 1 — (1— pc) 1 / 4 
for the center of a quarter sphere at a (d — 2)-dimensional boundary, and so on. It is only slightly harder to compute 
the covering probability for points in the bulk at distance p from a [d — l)-dimensional hypersurface, 



Pbulk(p) = 1 - cxp(-evj l (-r,p,r)/V d r d ) = 1 - (1 - p c )V$ {-r,P,r)/V d r* 



(24) 
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bulk covering boundary covering either covering 




FIG. 5: Probability of covering a point at a distance p from a (d— l)-dimensional boundary with a sphere from: left — a random 
bulk covering with pc = 95%; center — a random boundary covering with ps = 77%; right — either covering. As in the main 
text, here r is the covering radius. 

This function is shown in the left panel of Fig. [5] for pc = 95% and d = 2-20. 

There are two simple strategies to correct for boundary effects: extending the bulk covering across the boundary, 
and adding a supplementary covering on the boundary. We now examine each in turn. 

A. First strategy: Extend the bulk covering across the boundary 

Looking at the left panel of Fig. [5j we see that the covering probability for points in the bulk decreases significantly 
only at distance p < r/2 from the (d— l)-dimensional boundary. If we choose to extend the bulk covering to a distance 
r/2 beyond the boundary, we see that the total number of covering points must increase by a factor ~ (1 + (S/V) xr/2), 
where S is the hypersurface of the covered hypervolume V, at least for sufficiently smooth hypersurfaces. At any given 
dimension, the ratio S/V achieves the minimum possible value for a hypersphere, where S/V = d/R (for a hypercube, 
S/V — 2d/L). Thus the increase in the number of points incurred by extending the covering has a lower bound 
~ (1 + dr/D), where D is the linear spatial extent of the covered volume. The less symmetric the covered volume, 
the more important the boundary effects, which can be dominated by the dimension of shortest extent; of course if 
this extent is much smaller than the covering radius, the covering problem is effectively one of a lower dimension. 

We tested these theoretical predictions a la Monte Carlo, by repeatedly generating random coverings for the 
(i-dimensional hypercube. For all the tests in this paper we employed a Mersenne-twister pseudorandom number 
generator [43], which has extremely long period and is expected to generate uncorrelated sequences of points (i.e., 
d-uples of reals) up to 623 dimensions. (Linear congruential generators [25], by contrast, tend to generate d-uples 
that lie on hyperplanes, which would affect covering performance much as it does for periodic lattices.) We generated 
coverings of 

^ = 09r% d ° m (l + 2^(p^) (25) 

random points, uniformly distributed in the hypercube [—6,1 + S] d (with 5 — r or r/2); we then estimated the 
covering fraction by drawing 10,000 (target-signal) points uniformly across the hypercube [0, l] d . The resulting Ne 
and covering fractions are shown in Table [TV] For each d, we chose r so that a bulk-only covering would have between 
100 and 100,000 points (which are representative values for the template banks used in GW searches). Clearly the 
number of points in the extension across the boundary grows dramatically with d. 

This of course would be true also for periodic-lattice coverings, which have "overflow." Consider for instance the 
covering of [0, l] 8 with Z 8 and r — 0.3. Because 

(^W) * 453 (^059^) * 17 ' 0U = (3 - 38)8 ' (26) 

we need to choose between a covering of 3 8 = 6, 561 points with poor coverage near the sides, or a much more 
expensive covering of 4 8 = 65, 536 points. Boundary effects are certainly less important when r is smaller, but the 
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d 


Covering sphere 
r V d r d 


Bulk + border, 8 — r 
N cov.% 


Bulk + border, 8 = r/2 
N cov.% 


Bulk only 
N cov.% cov.% (theory) 


1 


10" 3 


2.00 x 10~ 3 


1,500 


95% 


1,499 


95% 


1,497 


95% 


95% 




0.01 


2.00 x 10~ 2 


152 


92-96% 


151 


92-96% 


149 92- 


-96% 


95% 


2 


0.01 


3.14 x 10~ 4 


9,920 


95% 


9,727 


95% 


9,535 


95% 


95% 




0.05 


7.85 x 10~ 3 


461 


~ 95% 


420 


~ 95% 


381 ~ 


94% 


94% 


3 


0.05 


5.24 x 10~ 4 


7,615 


95% 


6,623 


95% 


5,721 


93% 


94% 




0.10 


4.19 x 10~ 3 


1,235 


95% 


951 


95% 


715 


92% 


92% 


4 


0.05 


3.08 x 10" 5 


142,207 


95% 


118,062 


95% 


97,129 


93% 


93% 




0.10 


4.93 x 10~ 4 


12,588 


95% 


8,887 


95% 


6,070 


91% 


92% 




0.20 


7.90 x 10~ 3 


1,457 


95% 


786 


95% 


379 ~ 


88% 


90% 


6 


0.15 


5.89 x 10" 5 


245,650 


95% 


117,718 


95% 


50,892 


87% 


90% 




0.20 


3.31 x 10~ 4 


68,201 


95% 


27,046 


95% 


9,057 


84% 


89% 




0.30 


3.77 x 10" 3 


13,341 


95% 


3,838 


94% 


795 


78% 


87% 


8 


0.30 


2.66 x 10~ 4 


483,175 


95% 


91,768 


95% 


11,249 


74% 


86% 




0.40 


2.66 x 10~ 3 


124,112 


95% 


16,621 


95% 


1,126 


66% 


85% 




0.50 


1.59 x 10~ 2 


48,372 


95% 


4,842 


93% 


188 ~ 


57% 


84% 



TABLE IV: Number of additional templates and resulting covering fraction with the covering-extension strategy for a [0, l] d 
hypercube, as a function of dimension d and covering-sphere radius r. The required covering fraction is 95% for all runs. 
Columns 4 and 5 (6 and 7) show the number of covering points and the achieved covering fraction, to 1% accuracy, for a border 
of width 8 — r (8 = r/2). Columns 8 and 9 show the same information for a bulk-only covering; column 10 shows the theoretical 
covering fraction for a bulk-only covering obtained by considering only (d— l)-dimensional effects, which suggests that (d— 2)- 
and lower-dimensional boundaries become important as d increases. See the Appendix for details about numerical methods. 



overall numbers are much larger. For instance, for r — 0.05 we need to choose between 20 and 21 points to each side. 
Now, (21/20) 8 ~ 1.48, but 20 s = 2.56 x 10 10 . 

To apply the extension strategy on a curved manifold, we need to determine how far the covering must be extended 
in coordinate space to achieve the proper bordering, which requires knowledge of the metric. For GW template banks, 
we need also to worry that a waveform with parameters outside the original region of interest may not be physical, 
or may not even exist. For this reason algorithms such as the square-lattice placement of Ref. [17 place templates on 
the boundary first, and then in the bulk. This is just the second strategy that we consider next. 



B. Second strategy: Add a lower-dimensional covering on the boundary 

This scheme is generally easier to implement than the extension across boundaries, and it is a natural fit for the 
metricless, mesh-based placement methods described later in this paper. Now, how many points must be placed on a 
(d — l)-dimensional boundary to achieve uniform covering throughout the bulk? To answer this question, assume we 
lay down a boundary covering with covering probability ps, and compute the probability of covering a generic bulk 
point with at least one sphere from the boundary covering: 

PsuAp) = 1 - (1 - psf-WW- 1 "' for p < r; (27) 

this function is shown in the center panel of Fig. [5] for p$ — 0.78, and can be obtained easily by realizing that a 
(d — l)-dimensional boundary covering with radius r becomes a covering with radius r' on a (d — l)-dimensional 
surface parallel (at least locally) to the boundary. The new radius r' can be determined using the Pythagorean 

withV R = V (d ^ 1) (r'Y d - 1 \ 

Close to the boundary, the combined covering probability from the bulk and the boundary coverings will then be 

Pbulk+surf (p) = 1 - (l - Pbulk(p)) (l - Psurf (/»)) • (28) 

If we choose ps = 1 — exp(— 0/2), the points on the boundary get the same covering probability as those in the bulk, 
Pbuik+surf (0) = 1 — exp(— 0) = pc- Close to the boundary, 

Pbuik+suxf (p) = 1 - exp(-9 C{p)), with C(p) = V} (-r, p, r)/V d r d + i(l - (p/r) 2 )^/ 2 . (29) 



theorem (see Fig. 4b, and the relation between the covering probabilities using Eq. (15 
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As expected, Pbuik+surf is Pc for P = r i an d larger for < p < r. 

In principle it is possible to apply the same strategy to boundaries of dimension lower than d — 1, but the best 
arrangement will vary from case to case. Consider for example the 1-dimensional edges of a 3-dimensional cube: 
after placing a random covering in the bulk and on the 2-dimensional sides as we have just outlined, the covering 
probability on the edges is 1 — (1 — P3)(l — P2) 2 = 1 — cx p(~ (3/4)0), where P3 = 1 — cxp(— 0/4) is the probability of 
an edge point being covered by at least one point from the bulk, while p-2 = 1 — cxp(— (0/2)/2) is the probability of 
it being covered by a point from one of the sides shared by the edge. So an additional one-dimensional covering must 
be placed on the edge with pe = 1 — exp(— (0/4)) to achieve a covering probability 1 — exp(— 0) on the edges. Notice 
that to set p^we have used the information that two sides that share an edge are orthogonal, which will not be the 
case in general. To apply the boundary-covering strategy on a curved manifold, we need to compute the determinant 
of the metric on the hypersurface in terms of the d — 1 parameters that define it. 

How do the two strategies compare? For a [0, l] d hypercube, the points on the (d— 1) -dimensional boundary covering 
would be 



, 0\ 2dl d 



to be compared with the extension points 



_ l d {l + 2{r/2)/l) d ~l d ^ drl d ^ 

E ~ ^bulk - y dr d i 31 ) 

(where we assume r <C I). Thus Ne/Nq = V^-i/Vd, which is ~ 1 to very high d. This result is correct also for more 
general geometries, as long as the bordering volume scales as the hypersurface times r/2. However, the two strategies 
will in general have different distributions of r w (the radius of the largest empty sphere that can be fit among the 
covering points), since these depend on d. 



VI. SELF-AVOIDING COVERINGS 



As pointed out by Messenger and colleagues [25] and as discussed above, for sufficiently high d, random partial 
coverings become more efficient (i.e., have lower 0) than the best-known periodic lattices. In addition, they conform 
naturally to curved manifold geometries, given only a knowledge of the metric's determinant. However, each point in 
a random covering is placed independently of every other point, so it is natural to ask whether random coverings can 
be improved by making them self-avoiding, while preserving their random, unstructured character. In this section 
we consider two ways to introduce self-avoidance: the use of quasirandom sequences and the stochastic schemes that 
accept each random draw only after considering its distance from the points already accepted into the covering. 



A. Quasirandom sequences 

Quasirandom sequences, also known as low-discrepancy sequences 128,, are designed to cover multidimensional 
regions more uniformly than pseudorandom rt-tuples, although they may not appear as "random" (i.e., their deter- 
ministic nature is manifest). They are sometimes used to improve the convergence of multidimensional Monte Carlo 
integration. A conceptually simple example is Halton's sequence [28] . shown in Fig. [6] Quasirandom sequences can 
be used as slot-in replacements for pseudorandom numbers in random coverings. 

Effective thickness. Numerical experimentation with the more sophisticated and widely adopted Sobol sequences 
[31], shows that their effective thickness in (bulk) E d does significantly improve on the random-covering value, but it 
grows closer to the latter as d increases (see Table |v| and the left panel of Fig. [7|. Indeed, the quasirandom thickness 
may approach the random-covering value asymptotically as d — ¥ 00, but a proof (or disproof) remains to be found. 
In addition, the thickness of quasirandom coverings appears relatively stable under dilations along one axis. 

Curved manifolds and boundary effects. To use quasirandom coverings on curved manifolds, their density can 



IV 



be locally modulated in proportion to \J\gki{\)\, using a variant of the rejection method discussed in Sec 
Unfortunately, such a process tends to destroy the self-avoidance of quasirandom sequences, resulting in effective 
thicknesses close to the random-covering values (see right panel of Fig. [7]) . Thus the utility of quasirandom coverings 
appears limited to flat or almost flat manifolds. As for boundary effects, the geometrical analysis of Sec. [V] still 
applies, so the same covering strategies can be carried over, using empirically determined relations between thickness 
and covering percentage such as those plotted in the left panel of Fig. [7| 
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FIG. 6: Halton's quasirandom sequence. Left: to obtain the n-th number in the sequence, write n as a base-& number (here 2), 
reverse the digits, add a radix point at the left, and interpret the resulting string as a number € [0, 1). This algorithm fills a 
succession of finer and finer Cartesian grids, spreading out the points maximally on each, since the most rapidly changing &-ary 
digit of n controls the most significant digit of the placement [28]. Right: the combination of base-2 and -3 Halton sequences 
fills the unit square more uniformly than the same number of (pseudo-) randomly placed points. 



d 


©90% 


©95% 


©99% 


2 


1.6 


2.0 


2.9 


3 


1.7 


2.1 


3.0 


4 


1.8 


2.2 


3.2 


6 


1.9 


2.5 


3.6 


random 


2.3 


3.0 


4.6 



TABLE V: Numerically determined effective thickness of Sobol quasirandom coverings of bulk Euclidean space (no boundary 
effects). See Appendix for details on the numerical experiment. The same setup was used for Fig. IT] 



B. Stochastic placement 



Another approach to enforcing self-avoidance in random coverings is to assemble a covering by accepting or rejecting 
each random draw depending on already accepted points; such algorithms have been called stochastic. Harry and 
colleagues [26] describe a destructive process whereby an overdense covering is first generated; a point P is then 
chosen randomly from it, and all points closer to P than a distance are removed (with r the covering radius, and 



d 7 


©80% "80% 


©90% a 90% 


©95% C*95% 


©99% Q 99% 


2 0.5 
0.6 
0.7 
0.8 
0.9 
1.0 
1.1 


1.41 1.2 
1.33 1.3 
1.24 1.4 
1.20 1.6 
1.12 1.9 
1.07 2.2 
1.02 2.7 


1.90 1.3 
1.76 1.5 
1.61 1.7 
1.52 2.0 
1.42 2.4 
1.32 3.1 
1.24 4.4 


2.33 1.4 
2.13 1.6 
1.93 1.9 
1.77 2.4 
1.62 3.1 
1.50 4.5 
1.41 w 8 


3.10 1.7 
2.78 2.1 
2.46 2.6 
2.17 3.6 
1.97 5.8 
1.80 w 11 
1.66 w 60 


3 0.8 
1.0 


1.3 1.5 
1.2 2.2 


1.7 1.8 
1.5 3.2 


2.1 2.2 
1.8 4.6 


2.8 3.0 
2.2 w 12 


4 0.8 
1.0 


1.4 1.35 
1.3 2.2 


1.9 1.5 
1.7 3.1 


2.4 1.75 
2.0 4.5 


3.4 2.3 
2.7 w 12 


random 


1.61 1 


2.3 1 


3.0 1 


4.6 1 



TABLE VI: Effective thickness of stochastic coverings of bulk E d , for d = 2,3,4 and for different choices of covering fraction 
and acceptance radius 7. The factor a is the ratio between the number of proposed and accepted points. For comparison, the 
last row shows the effective thickness of random coverings, for which 7 and a are formally and 1. 



15 




FIG. 7: Left: effective thickness as a function of covering fraction for Sobol quasirandom coverings of bulk E d , for d = 2, 4 and 
6, as compared to the theoretical expectation for random coverings. (As a control, we used the same "experimental" setup to 
evaluate the random-covering thickness for d = 2, 4 and 6, which matches the theoretical expectation very accurately, and is not 
shown here.) Right: deterioration of quasirandom-covering thickness on curved manifolds. We estimate the effect of applying 
a rejection method in the presence of a varying \/\gki\ by extracting random subsets of 1/2, 1/5, and 1/10 of all points from 
a Sobol covering of E 2 , and evaluating the resulting thickness. Already for a range of variation ~ 10 in yjfffcTf, the thickness 
grows very close to its random-covering value for the same covering fraction. 



7 G [0, 2]); this is repeated until all remaining points have relative distances larger than jr. The resulting covering has 
slightly better thickness than a random covering with the same covering fraction. By contrast, Babak |24j implements 
a constructive process that begins drawing random templates, and rejects all new points that are closer than a given 
distance from already accepted points. In our analysis below we set this distance to the multiple j r of the covering 
radius. (An important difference between the methods of Refs. [35] and [H] is that Harry and colleagues use distances 
computed from a local metric, whereas Babak uses A distances, so his covering spheres can span multiple disconnected 
regions across parameter space.) 

Equivalence of destructive and constructive processes. The constructive and destructive processes are exactly 
equivalent: it is possible to formulate them in such a way that for the same input (a random sequence {P^'} 
of Nd candidate points) they would produce the same covering. Let us see why. In the destructive process, let 
|p0) } be the points in the overdense covering, and consider them for inclusion in the thinned covering in the order 
j = 1,2,..., Nd- In the constructive process, use the same order to evaluate each point for inclusion. In both cases, 
once a is accepted, its presence rules out accepting all P^ (with k > j) closer than jr to P^\ It is immaterial 
whether all such P^ are discarded immediately (as in the destructive process) or as they are "called up" (as in the 
constructive process). Thus the final set of accepted points is the same in both cases. Furthermore, the set of all the 
distances that must be computed is the same: it consists of the distances between each accepted P^ and all accepted 
PW> with k > j, and of the distances between each discarded and all accepted up to the first Pv ' closer 
than 7 r. 

Effective thickness. The choice of 7 is important. Setting 7 = 2 yields a solution of the sphere- packing problem 
|40j . so for 7 close to 2 not all covering fractions can be achieved, because there is a limit to the number of solid 
impenetrable spheres that can be fit in a given volume. A key quantity that can be determined empirically is the 
average number of draws a that are needed to accept a point: it is defined as a = Nd/N, wh ere Nd is the total number 



of draws, and N ~ ® s x%V m / (Vdr ) is the number of points of the final covering. Table VI shows the empirically 
estimated thickness and the a factor for stochastic coverings of bulk E 2-4 at different values of 7. The thickness 
improves (but a increases) for higher 7; the gains are diminishing at higher covering fractions, and at higher d. This 
confirms the general tendency that in higher dimensions it becomes progressively harder to improve the thickness of 
pure random distributions. 

Instead of keeping 7 constant, it is also possible to vary it while points are being added to the covering; in particular, 
by decreasing 7 oc N a (where N a is the number of points accepted so far) we can keep the number of attempts 
needed to accept a new point roughly constant. This can be understood as follows. The probability of a new point 
of being accepted is proportional to the fraction of the manifold left uncovered by the covering spheres of radius 7 r, 
and therefore to 1 — pc, where pc is the (average) covering fraction attained so far for a covering radius jr. Since 
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d a 


©80% 


©90% 


©95% 


©99% 


2 2.0 


1.15 


1.57 


1.99 


2.91 


3.0 


1.05 


1.42 


1.74 


2.50 


4.0 


1.01 


1.35 


1.66 


2.32 


10.0 


0.94 


1.19 


1.43 


1.89 


3 2.0 


1.23 


1.70 


2.14 


3.20 


3.0 


1.14 


1.55 


1.95 


2.9 


4.2 


1.11 


1.50 


1.87 


2.6 


12.0 


1.05 


1.39 


1.68 


2.3 


4 2.0 


1.3 


1.8 


2.3 


3.4 


4.5 


1.2 


1.7 


2.1 


3.0 


15.0 


1.2 


1.6 


2.0 


2.8 



TABLE VII: Effective thickness of scale-free (variable-7) stochastic coverings of bulk E d , for different choices of covering fraction 
and acceptance factor a. 



Nail r) d Vd — 0(pc)Vm, the covering fraction can be kept constant by changing 7 so that N a ^ d remains constant. 
The initial value of 7 is arbitrary, but it must correspond to an initial sphere volume (7 r) d Vd ~ V m if we are to 
produce a sequence with thickness Q™% appreciably different from a pure random covering. 

The thickness of such coverings is slightly worse than the thickness of the constant-7 coverings that yield the same 
average a (see Table VII). However, variable-7 coverings are interesting because they are scale free, just as random 
and quasirandom coverings: that is, they are independent of the covering radius used to build them. Thus they can 
be produced in advance and stored, and then used as a poor man's version of quasirandom sequences. To generate a 
covering with covering fraction X% and covering radius r, we would simply use the first N = <d™%V m /(Vdr d ) points 
in the sequence. However, their thickness is not stable under dilations, but converges to the thickness of pure random 
coverings. 

Computational cost. In principle we obtain the greatest covering efficiency by setting 7 as high as allowed by the 
desired covering fraction, and accepting a correspondingly large a. In the practice of matched filtering, however, we 
must balance the computational cost of placing a template with the cost of using it to filter the detector data. The 
latter is proportional to N, but the former is proportional to the number of distances that need to be computed in 
the stochastic process, which is < aN 2 /2. (This estimate includes N x N/2 distances between accepted points, and 
(a — 1)N x N/2 distances between discarded and accepted points. The second number is an upper limit because each 
discarded point could have been eliminated after comparing it with one of several accepted points; but the number 
of these neighbors is generally small compared to N.) The cost of computing distances may be negligible if these can 
be derived reliably from an analytic metric, or from a numerically obtained metric that is constant across H. On the 
other hand, the cost may be considerable if each distance requires the actual generation of templates, either to take 
their A distance directly, or to compute numerical metrics at different points in %. If the placement cost is dominated 
by generating templates as opposed to computing with them (e.g., taking their A distances), then its scaling can be 
attenuated to ~ aN by storing all templates as they are generated. Alternatively, the total number of distances that 
need to be computed can be made to scale as ~ d a N if it is possible to compare source parameters to decide which 
templates are likely to be neighbors. It is also true that placement cost is not an issue when detectors are stable 
enough that noise can be considered stationary [remember that noise levels affect distances through Eq. Q] ; in that 
case, the bank can be placed once and reused across a long dataset, so that the total computational cost would be 
dominated by the filtering. 

Curved manifolds. It would appear prima facie that the rejection process implicit in stochastic algorithm could 
replace the local density modulation of random draws needed by random coverings to cover curved H uniformly. This is 
correct, but there are two caveats. First, the stochastic process can become very inefficient on curved manifolds, since 
the average number of draws needed to accept a point is not a, but an average of a x (max^. \/\gki{K)\)/ \/\9ki(^i)\ 
across T-L . Second , the covering fraction of the final covering will not be uniform, but it will slightly favor the regions 
where \/|<?/d(Ai)| is lowest. Thus, if the determinant of the metric is available across H, we recommend combining 
the stochastic process with the nonuniform generation of random points. 
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signal Riemannian 




FIG. 8: Mapping of a parameter-space triangulation T onto the signal manifold H. If the A edge lengths are smaller than the 
characteristic curvature scale S of the manifold, the image of T will lie close to H. 

VII. METRICLESS, MESH-BASED RANDOM COVERINGS OF CURVED MANIFOLDS 

Throughout this paper, we have assumed that we could transfer the covering strategies formulated for E d to curved 
manifolds % using the pointwise knowledge of ^/Iffjfcl- (A possible exception are "straight" stochastic coverings of 
curved curved manifolds, which however can be very inefficient, as discussed above.) We now turn to the case where 
•v/f^fe] is not available. In the context of template placement, this may happen when we do not have a (semi-)analytical 
expression for the metric, because the waveform equations cannot be differentiated with respect to source parameters 
(e.g., the waveforms may be generated from the numerical solution of differential equations); the next recourse would 
be to take the numerical derivatives of Eq. ([3| for small source-parameter displacements, but this too can prove 
difficult because of numerical noise or computational cost. 

In this section we show how a discrete data structure consisting of a triangulation of coordinate (parameter) 
space and of the A distances measured along the triangulation's edges can be used in lieu of a metric to build 
properly density-modulated random coverings. (By contrast, Beauville and colleagues [5T] describe the use of a 
refined triangulation to interpolate equal-distance contours across two-dimensional parameter space to guide the 
placement of a locally hexagonal lattice.) 

Triangulations. To triangulate coordinate space, we decompose it into simplexes (the <i-dimensional analogues 
of triangles) such that their union covers the coordinate-space region of interest, and their intersection has zero d- 
dimensional volume. If the triangulation T is dense enough that the Riemannian distances measured in % along 
the triangulation edges are smaller than a scale 5 at which the curvature of % is negligible, then the image of T in 
M. Nhi & follows H closely (see Fig. |8|, and the Riemannian distances are approximated well by the A distances between 
neighboring vertices in T. In the following, we shall make this crucial assumption. We shall also assume that y]ftfe| 
is almost constant across each triangle; this condition depends of course on the choice of coordinates as well as the 
intrinsic geometry of %. 

By construction, we may now approximate the volume of a region of T-L by the sum of the volumes of all simplexes 
within the region, as computed using the Cayley-Menger determinant formula [35] and the A edge lengths. The 
triangular inequality Q for A guarantees that all simplexes have well-defined volumes. Furthermore, the ratio between 
the A volume of a simplex in H and its Euclidean volume in coordinate spaces approximates yf^fcj, although in 
practice this calculation can suffer from numerical noise. 

More formally, we see that the triangulation, together with the A edge lengths, carries the same information as the 
metric: if we assume that gjf. is constant across a simplex, consistently with our assumption of local flatness, we can 
recover gjk by solving the distance equations 

A A =g jk A£ A A£ A , A=l,...,d(d + l)/2, (32) 

where A enumerates the simplex edges, the A^ are their A lengths, and the Al A are the components of vectors that 
lie along the edges in coordinate space. 

Random coverings. Thus armed with an estimate of the global volume of the region of interest in H, and of the 
local \f\g~jk\, we are now able to generate a properly density-modulated random covering of H by a rejection method. 
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There is however a better algorithm that achieves the same result without discarding any random draw and without 



actually inverting Eq. (32 1 to compute the full metric: 



• To draw each new covering point, randomly select one of the simplexes in T in such a way that the probability 
of choosing each simplex is proportional to its A volume, then pick a point randomly in coordinate space within 
the chosen simplex. 

• To select the simplex, form an iVy-dimensional vector given by the cumulative sum of the A volumes of all 
simplexes (arbitrarily ordered), 



Wj = Y t Vj, (33) 



i=i 



then draw a random number x uniformly distributed in [0, Wjv T ], and choose the first simplex for which Wj > x. 

To pick a point uniformly within the chosen simplex, simple algorithms such as the following can be used. Begin 
by considering the d-dimensional unit simplex, which in K d+1 is embedded in the cZ-dimensional hyperplane 
Vi + 2/2 + • • • + Vd+i = 1, and has vertex coordinates y\ = b\ (for a, b = 1, . . . , d + 1). A random point in the 
unit simplex can be generated by drawing d+ I random numbers r a uniformly in [0, 1], and combining them as 
x a — ln(r a )/(ln(ri) + ■ • • + In(r<j+i)). This point can then be mapped to a point in an arbitrary d-dimensional 
simplex by the affine transformation Xj — xivj + • • • + Xd+iVj 1 , where the u" are the coordinates of the a-th 
vertex of that simplex. For details see Refs. [4"r?l 147] . 

Locally, this random covering is uniform: the assumption that the coordinates are distorted only by an affine trans- 
formation on the typical length-scale of a simplex edge ensures that \f\gjk\ is locally constant. Globally, by drawing 
points with density proportional to the A volume of simplexes, we correct for the variation of at larger scales. 

We can even estimate how many points we should draw to achieve a covering X% for covering radius r: from Eq. 
(|8j) that is of course N — 0™y^ om W jv t / '{VdT d ) . To correct for boundary effects, we may generate boundary coverings 
guided by the lower-dimensional triangulation consisting of the faces of the full-dimensional simplexes that lie on the 
boundary. 

Generating triangulations. It is of course pointless to discuss the advantages of triangulation-guided coverings if we 
cannot prescribe a convenient procedure to build appropriately dense triangulations of coordinate space. We propose 
a solution based on the incremental refinement of an initial sparse triangulation, which may be generated randomly. 
The refinement can be stopped when the triangulation satisfies a criterion based on the A edge lengths (for instance, 
if we have an estimate of the curvature scale S, we may require that all edge lengths be safely below it). The scale 
5 can be seen as a tuning parameter: the validity of an estimate for 5 can be checked a posteriori by measuring the 
covering fraction of the final coverings and its variation across parameter space. 

The actual details of the refinement process depend on what kind of triangulation we maintain. Consider for instance 
the Delaunay triangulation of a set of points |28j . This is the unique triangulation such that the circumsphere of each 
simplex contains no other point; it has the property of (maximally) avoiding "skinny" simplexes with small angles, 
and therefore it is a good choice to model terrains (for d — 2) or other hypersurfaces given a set of sample points. To 
refine a Delaunay triangulation, we can iteratively choose one of its simplexes on the basis of its A edge lengths, or 
of its A volume, place a new point at the barycenter, or randomly within the simplex, and retriangulate (see Fig. |9]). 
Efficient incremental algorithms exist that can adjust the triangulation to link the new point while preserving the 
Delaunay property (see, e.g., [28]). 

We have experimented with modifying these algorithms so that they would create triangulations that are Delaunay 
with respect to A-wise (rather than coordinatewise) circumsphcrcs. Such triangulations minimize the number of 
simplexes needed to approximate % faithfully but they suffer from a chicken-and-egg problem, because the local 
assessment of the Delaunay property is only reliable when the A edge lengths arc already below the curvature scale 
5, which was the whole point of refining the triangulation. An alternative to iterative Delaunay triangulations are 
longest-edge partition algorithms |48) . which refine an initial triangulation by iteratively placing a new point on the 
longest one-dimensional edge, and dividing all the simplexes that share that edge. Again it may be useful to evaluate 
the longest edge with respect to A distances. 

Refined-triangulation coverings. We note that it may be possible to use the very points of a refined triangulation 
as the points of a covering: in this case we would want to stop the refinement by comparing the population of edge 
lengths with the covering radius r. We investigate such algorithms in a separate paper [22], but we note here that 
the numbers of edges in a triangulation (and thus the number of A distances to compute) grows very rapidly with 
d, so triangulation-guided random coverings are still the better option if S 3> r and if coordinates can be found such 
that the variation of \f\g~jk\ across scales ~ 5 is small. If 5 ~ r, both triangulation-based approaches would still 
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FIG. 9: Three phases in the refined triangulation of E 2 , performed by placing new points at the barycenters of existing 
triangles. 

require fewer distance computations than "straight" stochastic placement, but the additional bookkeeping needed for 
the triangulation itself could become overwhelming, especially in higher dimensions. 

VIII. CONCLUSIONS AND FUTURE PROSPECTS 

Past and current searches for modeled GW sources have largely relied on filtering detector data with carefully 
distributed banks of theoretical signal templates. Furthermore, even bankless Monte Carlo searches (as envisaged for 
the space-based detector LISA [32) can benefit from the exhaustive a priori modeling of the posterior probability 
surface made possible by homogeneously distributed banks. The notion of a Riemannian metric in parameter space 
[TlTfT3"] allows placement methods based on periodic lattices [IT] [TH] , which however are limited in practice to simple 
signal models with very few source parameters. Thus, future searches could greatly benefit from more generic and 
robust placement methods that are suited to signal models with complex parameter dependencies and with moderate 
number of parameters. 

Template placement for generic signal families can be seen as an instance of the sphere-covering problem in Eu- 
clidean and Riemannian spaces. Working from this angle, Messenger and colleagues [35] examined the promise of 
random coverings, while Babak [24] and Harry et al. [26] studied stochastic coverings that combine random draws 
with the enforcement of a minimum distance between pairs of points. In this article we have developed a deeper un- 
derstanding of both kinds of coverings: specifically, we have derived analytically the variance of the covering fraction 
for random coverings, and examined the effects of boundaries; we have studied the self-avoiding coverings generated 
by quasirandom sequences; we have proved the equivalence of Harry et al.'s destructive stochastic coverings with 
Babak's constructive variant, and considered their effective thickness and computational cost; last, we have proposed 
a general technique to distribute coverings on curved signal manifolds using only the distances between the points of 
a parameter-space triangulation, removing the need for the Riemannian metric, which may be difficult to obtain. 

Overall, our study confirms that randomized (random and stochastic) coverings compare very favorably to lattice- 
based coverings, especially for higher-dimensional parameter spaces, where randomized coverings provide greater 
simplicity and flexibility with comparable thicknesses. Furthermore, unlike lattices, randomized coverings generalize 
straightforwardly from Euclidean to Riemannian signal manifolds; the required modulation of local density can be 
achieved by computing the determinant of the metric, but also by metricless methods such as "straight" stochastic 
algorithms that compare A distances, and by triangulation-based algorithms. 

Indeed, stochastic and triangulation approaches may be combined fruitfully: if the signal manifold has significant 
foldings (i.e., distinct parameter-space regions that correspond to very similar signals with small A distances), a 
sufficiently refined triangulation-based covering would separately populate each duplicate region, and (as shown by 
Babak [53]) a subsequent stochastic stage could recognize the foldings and generate a list of nonlocal mappings. Such 
a list would be a very useful input to Monte Carlo searches that need to jump between isolated peaks on the likelihood 
surface [32]. Furthermore, as pointed out by Babak [49], the initial stage of Monte Carlo searches (before the chains 
latch onto a candidate signal) can be seen as filtering by yet another flavor of random banks, so some of the methods 
and estimates developed in this paper and in Refs. |25[ I26j . as well as our discussion of boundary effects, could be 
useful in that context. 

Among the topics that we would like to flag for future investigation are the distance statistics of random and 
stochastic coverings, and especially the distribution of the maximum distance A max ; the possible mitigation of bound- 
ary effects with stochastic coverings, which naturally overpopulate the bulk regions near the boundaries; and the 
broad class of triangulation-based algorithms. It would also be interesting to investigate whether the interpolation of 
SNR across lattice-based template banks [T5l [38] [39] can be extended to the products of randomized placement. 
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As a final message, we wish to convey our belief, formed through the numerical experimentation carried out for 
this work, that the holy grail of a generally applicable template-placement algorithm is likely to remain unattainable: 
even general strategies such as random and stochastic coverings must be chosen, adapted, and carefully tuned for the 
specific search at hand. In every case, we first need to ask: which signals? what noise? what computational resources? 
The answer to these questions will guide the solution of what is arguably a problem of engineering applied to science. 
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Appendix A: Numerical procedures 



In this section we briefly describe the numerical procedures used to empirically determine covering fractions and 
effective thicknesses throughout this paper. 

Computing the thickness of partial periodic-lattice coverings. We work with the Voronoi cell of a lattice vertex, 
defined as the locus of points that are closest to that vertex than to any other (e.g., a hexagon centered on each vertex 
for the hexagonal lattice). We cover the Voronoi cell with a uniform distribution of points, and collect the vector Aj 
of their distances from the center. The X-th percentile of Aj is then the radius of the sphere covering that achieves 
a covering fraction of X%, and the corresponding thickness is 

e x% = V d (A x% ) d /V cclh (Al) 

where the Voronoi-cell volume V^ e ii is given by the determinant of the generator matrix |40j . This technique was used 
to generate the numbers of Table |H] and Figs, [l] and [3j 

Verifying the covering fraction of random coverings with boundary effects. As we do in the other two procedures 
described below, we begin by laying a very dense set of M target points uniformly distributed across the region of E'' 
to be covered, typically the hypercube [0, l] d (which of course has V m — 1)- We then place a random covering of N 
points throughout [0, l] d , or throughout the larger hypercube [—6, 1 + S] d , with S — r/2 or r, and verify what fraction 
of target points are covered (i.e., lie at a distance < r from the closest covering point). This technique was used to 
generate Table |IV| 

Computing the thickness of quasirandom coverings in the bulk. Again we lay a dense target set in [0, l) d ; we then 
place a covering of N points throughout a larger region that contains the hypercube. We compute the vector Aj of 
distances from the target points to the nearest point in the covering, and define 

O s: = N'V d (Ai) d /V m) (A2) 



where N' is the number of covering points that fall within [0, l] d . By definition, Eq. (A2| gives the thickness of a 
covering of V m with TV' points and covering radius A,-; by our very experiment, if we set Aj to its X-th percentile, 
such a covering achieves a covering fraction equal to X%. The X-th percentile of 8j is then the thickness of that 
covering. Boundary effects are avoided if the covering region is larger than [0, l] d by at least max.; Aj on every side. 
This technique was used to generate the numbers of Table [V] and Fig. [7] 

Computing the thickness of stochastic coverings in the bulk. We lay a dense target set on [0, l] d , pick r and 7, and 
place a stochastic covering with those parameters over a larger region that contains the hypercube. We then examine 
the covering points one by one, and keep a running tally Cj~ of the number of target points that have been covered 
by the first k covering points. The thickness for covering fraction X% is obtained from 

&x% = k'V d r d /V m , (A3) 

where we find the k such that is X% of M, and set k' equal to the number of covering points (among the first 



k) that fall within the hypercube. This technique was used to generate the numbers of Table VI (By contrast, the 



technique described in the paragraph above was used for the scale-free stochastic covering of Table VII ) 
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The last two techniques are equivalent, except that the first requires the a priori choice of the number of covering 
points, the second of the covering radius. 
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